import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import ConvexHull

class Point:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z

    def __repr__(self):
        return f"Point({self.x}, {self.y}, {self.z})"

    def __sub__(self, other):
        return Point(self.x - other.x, self.y - other.y, self.z - other.z)

    def __add__(self, other):
        return Point(self.x + other.x, self.y + other.y, self.z + other.z)

    def __mul__(self, scalar):
        return Point(self.x * scalar, self.y * scalar, self.z * scalar)

    def dot(self, other):
        return self.x * other.x + self.y * other.y + self.z * other.z

class Plane:
    def __init__(self, normal, point):
        self.normal = normal
        self.point = point

    def distance(self, p):
        vec = p - self.point
        return self.normal.dot(vec)

    def intersection_point(self, p1, p2):
        line_vec = p2 - p1
        denom = self.normal.dot(line_vec)
        if abs(denom) < 1e-6:
            return None  # 平行，无交点
        d = self.distance(p1)
        t = -d / denom
        if 0 <= t <= 1:
            return p1 + line_vec * t
        return None

class Polygon:
    def __init__(self, vertices):
        self.vertices = vertices

    def add_vertex(self, vertex):
        self.vertices.append(vertex)

    def get_vertices(self):
        return self.vertices

    def is_empty(self):
        return len(self.vertices) == 0

    def get_coordinates(self):
        return np.array([[v.x, v.y, v.z] for v in self.vertices])

def sutherland_hodgman(subject_polygon, clipping_planes):
    clipped_polygon = subject_polygon
    for plane in clipping_planes:
        new_polygon = Polygon([])
        vertices = clipped_polygon.get_vertices()
        prev_vertex = vertices[-1]
        for curr_vertex in vertices:
            prev_distance = plane.distance(prev_vertex)
            curr_distance = plane.distance(curr_vertex)
            if prev_distance <= 0 and curr_distance <= 0:
                new_polygon.add_vertex(curr_vertex)
            elif prev_distance > 0 and curr_distance <= 0:
                intersection = plane.intersection_point(prev_vertex, curr_vertex)
                if intersection:
                    new_polygon.add_vertex(intersection)
                    new_polygon.add_vertex(curr_vertex)
            elif prev_distance <= 0 and curr_distance > 0:
                intersection = plane.intersection_point(prev_vertex, curr_vertex)
                if intersection:
                    new_polygon.add_vertex(intersection)
            prev_vertex = curr_vertex
        clipped_polygon = new_polygon

        if clipped_polygon.is_empty():
            break

    return clipped_polygon

def plot_polyhedron(ax, polyhedron, color='b', alpha=0.6, edge_color='black'):
    vertices = polyhedron.get_coordinates()
    if vertices.size == 0:
        return

    hull = ConvexHull(vertices)
    ax.scatter(vertices[:, 0], vertices[:, 1], vertices[:, 2], color=edge_color, s=50, label='Vertices')  # 顶点加粗

    for simplex in hull.simplices:
        ax.plot_trisurf(vertices[:, 0], vertices[:, 1], vertices[:, 2], 
                        triangles=[simplex], color=color, alpha=alpha, edgecolor=edge_color)

def plot_plane(ax, point, normal, size=1.5, color='gray', alpha=0.3):
    """绘制有限大小的裁剪平面"""
    d = -np.dot([normal.x, normal.y, normal.z], [point.x, point.y, point.z])

    # 限定平面的大小
    xx, yy = np.meshgrid(np.linspace(0, size, 10), np.linspace(0, size, 10))
    zz = (-normal.x * xx - normal.y * yy - d) / normal.z

    # 绘制有限平面
    ax.plot_surface(xx, yy, zz, color=color, alpha=alpha)

if __name__ == '__main__':
    # 定义一个立方体
    cube_vertices = [
        Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
        Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
    ]
    subject_polygon = Polygon(cube_vertices)

    # 定义一个斜裁剪平面
    plane_normal = Point(1, 1, 1)  # 平面方向为 (1, 1, 1)
    plane_point = Point(0.5, 0.5, 0.5)  # 平面经过 (0.5, 0.5, 0.5)
    planes = [Plane(plane_normal, plane_point)]

    # 执行裁剪
    clipped_polygon = sutherland_hodgman(subject_polygon, planes)

    # 打印裁剪结果的顶点
    print("裁剪后的顶点坐标：")
    for vertex in clipped_polygon.get_vertices():
        print(vertex)

    # 绘制原始和裁剪后的多面体
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # 绘制原始多面体
    plot_polyhedron(ax, subject_polygon, color='#fcbad3', alpha=0.6, edge_color='black')

    # 绘制裁剪平面
    plot_plane(ax, plane_point, plane_normal, size=1, color='#a8d8ea', alpha=0.5)

    # 如果裁剪结果非空，则绘制裁剪后的多面体
    if not clipped_polygon.is_empty():
        plot_polyhedron(ax, clipped_polygon, color='#3f72af', alpha=0.7, edge_color='red')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('3D Sutherland-Hodgman Clipping with a Plane')
    ax.legend()
    plt.show()
